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ABSTRACT 

A  Fourier  method  is  used  to  model  mountain  waves  that  have  nearby  turning  points  in  a  wind  jet.  In 
Fourier  space,  the  propagation  equations  are  solved  by  ray  theory.  To  correct  for  the  ray  singularity  at  a 
turning  point  without  time-consuming  special-function  evaluations,  the  ray  solution  is  linearly  interpolated 
across  the  breakdown  region.  The  Fourier  solutions  for  the  spatial  wavefield  are  compared  with  mesoscale 
model  simulations  in  two  cases:  two-dimensional  flow  over  idealized  topography  with  uniform  stratification 
and  a  sech-squared  wind  profile  and  three-dimensional  flow  over  the  island  of  Jan  Mayen  with  stratification 
and  wind  profiles  taken  from  radiosonde  measurements.  The  latter  case  reveals  the  partial  transmission  of 
trapped  mountain  waves  into  the  stratosphere. 


1.  Introduction 

Mountain  waves  propagate  to  great  heights  in  the  at¬ 
mosphere  (see  the  review  by  Fritts  and  Alexander  2003), 
at  times  reaching  the  mesosphere  and  lower  thermo¬ 
sphere  (Bacmeister  1993;  Eckermann  et  al.  2007).  The 
propagation  can  be  interrupted  by  wind  jets,  which  pro¬ 
duce  evanescent  layers  that  partially  reflect  and  partially 
transmit  the  waves  (e.g.,  Nault  and  Sutherland  2008, 
hereafter  NaSu).  Above  the  wind  jet,  the  partially  trans¬ 
mitted  waves  can  continue  to  grow  with  height  in  re¬ 
sponse  to  the  decreasing  air  density  and  can  make  an 
important  contribution  to  the  wavefield  at  higher  alti¬ 
tudes  (see  the  numerical  simulations  in  Eckermann  et  al. 
2006,  hereafter  EB06). 

Flere  we  use  a  Fourier  method  to  model  mountain 
waves  in  the  presence  of  a  wind  jet.  In  Fourier  space  the 


Corresponding  author  address:  Stephen  Eckermann,  Space  Sci¬ 
ence  Division,  Naval  Research  Laboratory,  Code  7646,  Wash¬ 
ington,  DC  20375. 

E-mail:  stephen.eckermann@nrl.navy.mil 
DOI:  10.1175/2008JAS2786.1 


theory  is  one-dimensional,  and  the  solution  can  be 
computed  efficiently  by  ray  theory,  although  ray  theory 
breaks  down  near  turning  points.  For  a  single  isolated 
turning  point,  the  solution  is  efficiently  calculated  with 
an  Airy  function,  but  in  the  cases  we  consider  here  there 
are  as  many  as  four  nearby  turning  points  in  the  wind 
jet.  The  solution  for  just  two  nearby  turning  points  in¬ 
volves  Weber  (or  parabolic  cylinder)  functions  and  is 
very  slow  computationally. 

We  aim  for  computation  times  of  minutes  for  the 
entire  Fourier  spectrum,  as  has  been  achieved  (without 
the  treatment  of  nearby  turning  points)  in  previous 
applications  of  Fourier  methods  to  mountain  waves 
(Broutman  et  al.  2003;  EB06).  Efficiency  is  important 
because  in  realistic  three-dimensional  cases  there  can  be 
many  thousands  of  Fourier  components  with  nearby 
turning  points  in  the  wind  jet. 

We  find  that  an  efficient  but  approximate  solution  in 
the  region  of  nearby  turning  points  is  obtained  by  lin¬ 
early  interpolating  the  ray  solution  across  the  turning- 
point  region.  No  special  functions  are  evaluated,  but  it  is 
necessary  to  have  a  reliable  estimate  for  where  ray 
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theory  breaks  down  near  the  turning  points.  This  indi¬ 
cates  where  to  start  and  stop  the  interpolation.  The 
criterion  is  discussed  in  section  3d. 

We  first  test  the  interpolated  ray  solution  for  selected 
Fourier  components  using  a  sech-squared  mean  wind 
profile  and  constant  stratification.  We  consider  cases  in 
which  the  curvature  of  the  wind  profile  is  ignored  and 
included  in  the  theory.  When  ignored,  the  maximum 
number  of  turning  points  per  Fourier  component  is  two. 
When  included,  the  maximum  number  is  four.  We  do 
not  have  a  four-turning-point  theory,  so  we  attempt, 
with  some  success,  to  model  four  turning  points  with  the 
two-turning-point  theory. 

We  then  calculate  spatial  solutions  for  the  mountain 
waves  by  a  Fourier  superposition  of  the  ray  solutions.  We 
start  with  two-dimensional  mountain  waves  generated  by 
idealized  topography,  a  sech-squared  wind  profile,  and 
uniform  stratification.  We  then  consider  three-dimen¬ 
sional  mountain  waves  generated  by  the  island  of  Jan 
Mayen,  with  stratification  and  wind  profiles  taken  from 
radiosonde  measurements  over  the  island. 

The  Jan  Mayen  case  has  been  examined  by  EB06.  Their 
Fourier  solution  did  not  include  the  treatment  of  nearby 
turning  points  and  did  not  accurately  reproduce  the  results 
of  a  mesoscale  model  simulation  at  heights  above  the  wind 
jet,  centered  near  10-km  altitude.  From  estimates  of  the 
integral  given  here  by  Q  in  (7),  these  authors  suggested 
that  there  would  be  significant  wave  transmission  through 
the  jet  region.  We  repeat  their  Fourier  calculation,  but  we 
include  the  treatment  of  nearby  turning  points.  We  find  an 
improved  comparison  between  the  solutions  of  the  Fou¬ 
rier  method  and  the  mesoscale  model. 

In  quantum  mechanics,  the  evanescent  region  be¬ 
tween  turning  points  is  called  a  potential  barrier,  and 
wave  transmission  through  a  potential  barrier  is  called 
tunneling.  Our  four-turning-point  case  is  an  example  of 
a  double  potential  barrier.  Analogies  to  quantum  me¬ 
chanical  tunneling  for  gravity  waves  are  discussed  in 
Sutherland  and  Yewchuk  (2004)  and  Brown  and  Su¬ 
therland  (2007).  These  authors  derive  analytic  wave 
solutions  using  piecewise  constant  profiles  for  the 
background  stratification  and  shear.  Further  work  on 
gravity  wave  tunneling  by  Nault  and  Sutherland  (2007) 
and  NaSu  is  based  on  the  numerical  integration  of  an 
equation  similar  to  (1)  below  and  will  be  discussed  in  the 
concluding  section  of  our  paper. 

We  formulate  the  problem  in  section  2,  derive  two- 
turning-point  solutions  in  section  3,  consider  the  effects 
of  wind  curvature  in  section  4,  and  introduce  other  ef¬ 
fects  needed  for  mountain-wave  simulations  in  section 
5.  An  idealized  case  with  a  sech-squared  wind  profile  is 
given  in  section  6.  A  realistic  case  for  Jan  Mayen  is  given 
in  section  7,  followed  by  a  discussion  in  section  8. 


2.  Formulation 

We  consider  stationary  mountain  waves  propagating 
in  a  background  that  depends  on  height  z.  The  mean 
buoyancy  frequency  is  N(z).  The  mean  wind  has  com¬ 
ponents  U{z),  V{z)  in  the  horizontal  directions  x,  y, 
respectively.  Initially,  the  Boussinesq  approximation  is 
made.  In  section  5c  we  add  the  usual  anelastic  density 
scaling  for  the  wave  amplitudes. 

The  mountain  waves  have  horizontal  wavenumber 
components  k,  l  in  x,  y  respectively,  and  intrinsic  fre¬ 
quency  co  =  —  kU  —  IV.  The  vertical  velocity  for  a 
particular  Fourier  component,  w(k,l,z),  satisfies 

wzz+m2(z)w  =  0,  (1) 

where  (see  Teixeira  et  al.  2004) 

m 2  =  k l(N2/u2  -  l)+a(kUzz+lVzz)/w  (2) 

and  k/,  =  (k2  +  l2)112.  In  the  anelastic  approximation,  a 
term  —1/4 H2d  is  added  to  the  right-hand  side  of  (2), 
where  Hd  is  the  density  scale  height  (e.g.,  NaSu).  We 
ignore  this  term,  which  is  small  for  the  horizontal 
wavelengths  of  importance  here  (10-50  km). 

The  value  of  a  in  (2)  will  be  set  to  unity  or  zero, 
respectively,  to  include  or  omit  wind  curvature  terms. 
When  a  =  0,  (2)  is  the  dispersion  relation  for  internal 
gravity  waves  and  m  is  the  vertical  wavenumber.  It  can 
be  shown  that  the  vertical  group  velocity  is  positive 
when  m  and  d>  have  the  opposite  sign,  and  thus  when  m 
and  kU  +  IV  have  the  same  sign.  Our  notation  is  that  m 
is  positive  (negative)  for  upgoing  (downgoing)  wave 
groups.  This  is  opposite  to  the  sign  notation  often  used 
for  gravity  waves,  but  we  have  chosen  it  to  simplify  the 
application  of  two-turning-point  theories.  These  were 
developed  in  other  contexts  and  are  typically  formu¬ 
lated  so  that  the  wavenumber  is  positive  for  waves  that 
are  incident  upon  and  transmitted  through  the  turning- 
point  region. 

The  topography  has  elevation  h  (x,  y).  The  linearized 
lower  boundary  condition  for  mountain  waves  is  [cf.  Eq. 
(5.2.4)  in  Baines  1995] 

w  =  —  icoh,  (3) 

evaluated  at  z  =  0,  where  h  (k,  l)  is  the  Fourier  trans¬ 
form  of  h  (x,  y). 

An  inverse  Fourier  transform  of  w(k,l,z)  gives  the 
spatial  solution: 


w(x,y,z) 


w(k,l,z)eL{kx+ly)dkdl. 


—oo 


(4) 
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Fig.  1.  Wind  profiles  used  in  the  present  paper:  (left)  the  sech-squared  profile  in  (26);  (right)  the  westerly  com¬ 
ponent  of  the  wind  obtained  from  radiosonde  measurements  and  used  (along  with  a  weaker  north-south  component) 
in  EB06.  Turning  points  for  a  particular  Fourier  component  are  indicated  at  heights  z\  and  Zz-  They  separate 
propagating  regions  (m2  >  0)  from  evanescent  regions  (m2  <  0). 


3.  Two-turning-point  theory 

We  first  consider  the  solution  for  Fourier  components 
with  two  turning  points.  Figure  1  illustrates  the  problem 
for  the  two  wind  profiles  used  in  this  study:  a  sech- 
squared  profile  and  the  wind  profile  measured  above 
Jan  Mayen.  The  turning  point  heights  z.\  and  zi,  with 
Zi  <  z.2-  are  functions  of  k,  /  as  well  as  of  U,  V,  and  N. 

We  examine  three  solutions  to  (1):  the  ray  solution,  a 
uniformly  valid  solution,  and  a  solution  by  numerical 
integration. 

a.  The  ray  solution 

The  ray  solution  for  (1)  is  derived  by  assuming  that 
the  vertical  wavenumber  m  is  a  slowly  varying  function 
of  height.  The  derivation  also  yields  a  solution  in  eva¬ 
nescent  regions,  where  m  is  imaginary,  but  il  diverges  at 
turning  points,  where  m  =  0.  We  will  refer  to  this  as  the 
ray  solution  even  though  it  includes  regions  of  evanes¬ 
cence.  We  might  have  called  it  a  phase  integral  solution 
(as  in  Froman  and  Froman  2002  and  Fleading  1962),  but 
this  term  is  less  familiar  and  does  not  reflect  our  use  of 
ray  theory  to  incorporate  wave  transience  (see  section 
5).  The  term  “WKBJ  solution”  is  sometimes  used,  but 
this  is  often  meant  to  imply  that  the  ray  solution  has 
been  matched  to  a  special-function  solution  near  the 
turning  point.  We  have  not  done  the  matching  but  will 
instead  compute  a  uniformly  valid  expression  from  the 
ray  solution  (see  section  3b). 

The  ray  solution  used  here  is  derived  in  Froman  and 
Froman  (1965,  1970).  We  define  the  phase  integrals 


<72  (z)  =  /  mdz ,  and  (6) 

JZ2 

Q=  f  \m\dz.  (7) 

Jzi 

Note  that  below  the  lower  turning  point,  z  <  Zi,  we  have 
from  (5)  negative  q1  but  positive  dq^ldz.  Thus,  in  the 
following  ray  solutions  we  associate  ±qi  with,  respec¬ 
tively,  upgoing  and  downgoing  wave  groups  below  the 
lower  turning  point. 


The  ray  solutions  are 

w,  =  Cm_1'/2e^?2+7r/4), 

(8) 

wr  =  C\m\~V2eQe-i^+^\ 

and 

(9) 

wi  =  Cm-1/2{e2Q+V)V2ei(q'+ 

77/4  +  20") 

5 

(10) 

where  C  is  a  constant  determined  by  a  boundary  condi¬ 
tion.  The  subscripts  t ,  r,  and  i  refer,  respectively,  to  the 
transmitted  wave  (z  >  z 2)  and  the  reflected  and  incident 
waves  (z  <  Z\).  The  phase  shift  2<r  can  be  ignored  except 
when  the  two  turning  points  are  close  to  each  other,  ft  is 
at  most  about  10°,  as  noted  in  Berry  and  Mount  (1972). 
An  expression  for  <x  is  given  by  Froman  and  Dammert 
[1970,  Eq.  (57)]: 

cr  =  0.5[<2*  log  Q*  —  Q*+  argr(0.5  —  iQ*)\,  (11) 

where  Q*  =  Qhr  and  I  is  the  gamma  function. 

Thus, 

w(z)  =  Wi+wr  for  Z<Z1,  and  (12) 
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w(z)=w,  for  z>z2.  (13) 

The  upper  boundary  condition  is  a  radiation  condition, 
satisfied  by  the  form  chosen  for  wt,  which  represents  an 
outgoing  wave. 

Note  that  the  ray  solution  for  the  transmitted  wave 
can  be  derived  without  matching  the  ray  solution  to  a 
special-function  solution  valid  in  the  turning-point  re¬ 
gion.  The  relation  between  the  ray  solution  above  the 
turning-point  region  (z  >  z2)  and  the  ray  solution  below 
the  turning-point  region  (z  <  Zi)  can  instead  be  found 
by  following  the  ray  solution  through  the  complex  z 
plane  on  a  path  that  avoids  the  turning  points.  This  use 
of  the  complex  plane  was  suggested  soon  after  the 
original  WKB  papers  were  published  in  the  1920s.  For 
more  history  and  explanation,  see  Froman  and  Frbman 
(2002)  and  Fleading  (1962). 

Between  the  two  turning  points,  where  m  is  imagi¬ 
nary,  we  use  the  evanescent  WKB  solution  [see  Froman 
and  Froman  1965,  Eq.  (9.20b)]: 

w(z)  =  C\m\~x^e~^  for  z\<z<z2-  (14) 

This  single-exponential  form  is  only  valid  in  the  limit  of 
well-separated  turning  points  and  represents  the  lowest- 
order  WKB  solution.  Like  the  ray  solutions  in  (8)-(10), 
solution  (14)  diverges  at  each  turning  point  and  is  in¬ 
accurate  at  distances  close  to  each  turning  point. 

b.  The  uniformly  valid  solution 

A  uniformly  valid  solution  can  be  derived  that  ap¬ 
proximates  the  ray  solution  far  from  the  turning  points 
and  is  also  accurate  near  the  turning  points.  For  a  single 
turning  point,  the  uniform  solution  involves  the  Airy  func¬ 
tion  (e.g.,  Berry  and  Mount  1972).  For  two  turning  points, 
the  uniform  solution  involves  the  parabolic  cylinder  func¬ 
tion  (Berry  and  Mount  1972;  Razavy  2003),  which  is  a 
solution  to  the  parabolic  cylinder  equation  (Abramowitz 
and  Stegun  1964,  hereafter  AS64,  section  19.16) 

Wz.z,  +  {z'2/A-a)W  =  Q.  (15) 

Flere,  z'  is  a  nondimensional  height  coordinate,  and  the 
two  turning  points  are  located  at  z!  =  ±2 a1!2.  Equation 
(15)  has  a  coefficient  that  is  symmetric  and  parabolic  in 
z' ,  the  simplest  form  of  our  original  Eq.  (1)  with  a 
smooth  coefficient  and  two  turning  points.  Although 
parabolic,  the  solutions  to  (15)  can  be  scaled  to  ap¬ 
proximate  the  solution  for  more  general  two-turning- 
point  coefficients  that  are  not  parabolic,  except  in  the 
close  vicinity  of  the  turning  points. 

Two  independent  solutions  of  (15)  are  given  by  AS64 
as  E  (a,  z')  and  its  complex  conjugate  E*  (a,  z')-  Only 


E  ( a ,  z')  satisfies  the  radiation  condition  at  z'  3>  a.  In 
this  limit,  we  have  [from  AS64,  Eq.  (19.21.1)] 

E{a,  z')~(2 /z')  1/2e;(2'2/4_fl  log (16) 

where  f>2  =  arg  T  (1/2  +  ia).  The  derivative  with  respect 
to  z'  of  the  imaginary  part  of  the  exponent  in  (16)  gives 
a  wavenumber  m  ~  z' 12  —  alz' .  Thus,  m  >  0  for  z'  3>  a, 
and  E  represents  an  upward-propagating  wave  group 
above  the  second  turning  point. 

For  negative  z' ,  we  evaluate  E  from  AS64’s  Eq. 
(19.18.3): 

iE(a,-z')=  -(1  +  e2m)U2E*{a,z')+e7raE(a,z').  (17) 

This  expression  is  for  z'  >  0,  so  E  on  the  left-hand  side 
above  has  a  negative  second  argument.  For  (15),  the 
integral  in  (7)  is  Q  =  i to.  Comparing  the  exponential 
terms  in  (17)  to  the  Q  terms  in  the  ray  solutions  (9)- 
(10),  we  identify  the  first  term  in  (17)  as  the  incident 
wave  and  the  second  term  as  the  reflected  wave. 

The  idea  of  the  uniform  solution  (e.g.,  Berry  and 
Mount  1972)  is  to  adjust  the  size  and  shape  of  E  (a,  z') 
so  that  it  approximates  the  solution  of  a  more  general 
two-turning-point  equation.  This  adjustment  involves 
only  quantities  that  are  used  in  the  ray  solution  (Berry 
1969),  such  as  the  phase  integral  f  mdz  that  appears  in 
(8)-(10). 

To  simplify  the  expressions  below,  we  define 

£  =  z!  /  2} '2  and  b  =  (2a)1/2.  (18) 

This  transforms  the  parabolic  cylinder  Eq.  (15)  to 

W(C  +  (£2  -  b2)W  =  0,  (19) 

which  Kravtsov  and  Orlov  (1999)  call  Weber’s  equation. 
The  turning  points  are  at  £  =  ±b. 

The  uniform  solution  for  (1)  has  the  form  (Kravtsov 
and  Orlov  1999) 

w(z)=DA(z)E[b,£{z)],  (20) 

where  D  is  a  constant  used  to  satisfy  the  lower  boundary 
condition,  and 

Mz)  =  [(<T2(z)  -  b2)/m2(z)}1/4,  (21) 

b  =  2Q/ir.  (22) 

The  mapping  £  (z)  is  defined  by  the  following  inte¬ 
grals.  [Note  that  this  is  the  mapping  £  (z),  where  z  is 
our  height  coordinate,  not  £  (z')  of  (18),  where  z'  is 
the  coordinate  in  the  parabolic  cylinder  Eq.  (15).]  For 

z  <  Zi, 
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J  b  [if2  -  b2)1/2d(  =  mdz. 


For  zi  <  z  <  zi. 


For  z  >  Z2, 


[  (b2  —  g2)1/2 d£  =  [  \m\dz. 
J-b  J  zi 


f  (£2  -  b2)l!2d£=  I  mdz. 
Jb  J  Zi 


(23) 


(24) 


(25) 


The  turning  points  z  =  Zi  and  z  =  Zi  are  mapped  to 
£  =  —b  and  £  =  b,  respectively.  The  integrals  on  the 
left-hand  side  of  the  above  three  equations  can  be 
evaluated  analytically,  but  the  integrals  on  the  right- 
hand  side  generally  need  numerical  solution.  Also, 
the  explicit  relation  for  f  (z)  needs  to  be  determined 
iteratively. 


c.  Solutions  for  a  sech-squared  wind  profile 

We  consider  the  wind  profile  shown  in  the  left  panel 
of  Fig.  1: 


U  =  t/o  +  t/i  sech2[(z  -  Zmf/L2},  (26) 

with  Uo  =  10  m  s_1,  U]  =  30  m  s_1,  z„,  =  10  km,  and  L  = 
3  km.  The  other  background  profiles  are  V  =  0  and  N  = 
0.015  s~ 1 . 

For  now,  we  ignore  the  wind  curvature  terms,  setting 
a  =  0  in  (1).  There  are  then  at  most  two  turning  points, 
symmetrically  located  below  and  above  z  =  zm. 

All  numerical  results  in  the  present  paper  were  pro¬ 
duced  with  Matlab  version  7.4.  The  phase  integrals  in 
the  ray  solution,  (5)-(7),  and  the  integrals  in  the  uni¬ 
form  solution,  (23)-(25),  were  approximated  by  the 
trapezoidal  rule.  The  iterative  solution  to  (23)-(25)  was 
obtained  using  Matlab’s  root-finding  routine  fzero.  The 
parabolic  cylinder  function  in  (20)  was  evaluated  using 
Matlab’s  generalized  hypergeometric  function  routine. 
The  numerical  solution  to  (1)  was  obtained  with  Mat¬ 
lab’s  ordinary  differential  equation  solver  ode23,  with 
default  error  tolerances  of  10  1  for  the  relative  error 
and  1CT6  for  the  absolute  error. 

Figure  2  shows  the  ray  and  uniform  solutions  of  (1), 
along  with  a  solution  by  numerical  integration  (see  be¬ 
low).  The  left  and  right  panels  are  for  a  horizontal 
wavelength  of  16.3  and  14  km,  respectively.  The  solu¬ 
tions  were  computed  at  100  evenly  spaced  heights  from 
0  to  20  km.  The  ray  solution  diverges  at  both  turning 
points,  whose  heights  are  denoted  in  the  figure  by  hor¬ 
izontal  dotted  lines.  The  computed  values  are  finite 


because  the  discrete  grid  points  miss  the  turning  points 
slightly. 

The  numerical  solution  to  (1)  was  calculated  with  an 
initial  condition  given  by  the  transmitted  ray  solution 
(8)  at  z  =  20  km.  This  was  then  integrated  downward  to 
z  =  0.  The  numerical  solution  separates  slightly  from 
the  uniform  solution  as  the  integration  proceeds 
downward  through  the  evanescent  region  between  the 
turning  points,  although  the  agreement  is  still  quite 
good  below  the  wind  jet.  We  repeated  the  numerical 
integration  with  much  smaller  error  tolerances,  but  the 
result  was  the  same.  The  small  discrepancy  between 
the  numerical  and  ray  solutions  may  thus  be  due  to 
errors  in  the  ray  approximation.  Higher-order  correc¬ 
tions  to  the  ray  solution  can  be  computed  by  the 
method  of  Froman  and  Froman  (2002)  but  were  not 
attempted  here. 

In  our  calculations,  the  ray  solution  was  by  far  the 
fastest  to  compute.  The  numerical  solution  took  over  10 
times  longer,  and  the  uniform  solution  took  over  100 
times  longer  (about  a  couple  of  seconds  on  a  1-GHz 
processor).  Most  of  the  computation  time  for  the  uni¬ 
form  solution  was  spent  evaluating  the  parabolic  cylin¬ 
der  function  E  in  (20),  although  the  iterative  solution  for 
the  mappings  (23)-(25)  by  itself  took  10  times  longer 
than  the  ray  calculation.  We  experimented  somewhat 
with  parameter  settings  in  the  Matlab  routines,  and  the 
above  comments  refer  to  our  best  times.  There  may  be 
much  faster  algorithms  for  evaluating  the  parabolic 
cylinder  function. 

In  the  Jan  Mayen  example  reported  below,  we  need 
to  compute  two-turning-point  solutions  for  thousands  of 
horizontal  wavenumbers.  To  do  this  with  reasonable 
speed,  we  use  the  ray  method  and  introduce  a  simple 
modification  to  correct  its  breakdown  near  the  turning 
points,  as  we  now  discuss. 


d.  The  interpolated  ray  solution 

We  first  identify  the  region  of  ray-theory  breakdown 
using  the  following  quantity,  which  measures  the  slow 
variation  of  the  vertical  wavenumber  m: 


1 

3 

( 1  dm\ 2 

1  d2m 

1  |2 
m 

4 

\m  dz  J 

2m  dz2 

1 

16/7I6 


5 


2 

—4m2 


d2m2 

dz2 


(27) 


(28) 


Ray  methods  are  asymptotically  valid;  that  is,  their  ac¬ 
curacy  improves  as  e  — >  0.  At  a  turning  point,  e  diverges 
because  neighboring  rays  with  slightly  different  values 
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Fig.  2.  Solutions  for  w  for  the  sech-squared  wind  profile  (26)  for  (left)  horizontal  wavelength  16.3  km,  zi  =  9.4  km, 
and  Z2  —  10.6  km  and  (right)  horizontal  wavelength  14  km,  zi  =  8.2  km,  and  z2  =  11.8  km.  Horizontal  dotted  lines 
indicate  the  turning-point  heights  z  i ,  z2.The  solutions  are  normalized  so  that  w  =  - 1  at  z  =  20  km,  and  the  real  part  of 
w  is  plotted.  Note  that  these  are  Boussinesq  solutions.  Anelastic  scalings  are  added  in  section  5c. 


of  m  cross  each  other.  This  results  in  the  divergence  of 
the  derivatives  in  the  above  expressions. 

The  above  expressions  are  taken  from  Froman  and 
Froman  [1965,  Eq.  (3.5)].  The  first  form,  (27),  is  also 
given  by  Einaudi  and  Flines  (1970).  The  second  form. 
(28),  is  in  terms  of  m2,  which  is  real  for  all  real  z. 

We  consider  ray  theory  to  be  sufficiently  accurate 
wherever  e  <  1,  and  we  leave  the  ray  solution  as  it  is  at 
those  heights.  At  heights  where  e  >  1,  we  replace  the  ray 
solution  with  linear  interpolated  values  across  those 
heights.  Let  z  =  r  and  z  =  s  be  heights  on  either  side  of  a 
region  where  e  >  1.  Our  modification  to  the  ray  solution 
is  then  given  by  the  linear  interpolation 

i (29) 

We  show  examples  of  the  interpolated  ray  solution 
(indicated  by  thick  dashed  lines)  in  Fig.  3.  The  two 
upper  panels  are  the  same  cases  as  the  respective  panels 
of  Fig.  2.  We  approximated  e  using  centered  differences 
for  the  derivatives  of  m2  in  (28). 

The  lower  left  panel  of  Fig.  3  shows  the  least  accurate 
result  we  found  using  the  interpolated  ray  method 
for  the  present  wind  profile,  ignoring  wind  curvature. 
Here  e  reaches  a  local  minimum  of  about  0.8  near  the 
middle  of  the  wind  jet  at  z  =  10  km  (bottom  right 
panel).  We  experimented  with  other  values  of  e  for 
determining  ray  theory  validity,  but  the  value  of  unity 
usually  led  to  the  best  approximation.  A  possible  im¬ 
provement  to  our  procedure  might  be  to  interpolate 
across  the  entire  evanescent  region  when  the  turning 
points  are  sufficiently  close. 


4.  Wind  curvature 

Next  we  include  the  previously  omitted  wind  curva¬ 
ture  terms  by  setting  a  =  1  in  (2).  The  difficulty  is  that 
some  Fourier  components  now  have  four  nearby  turn¬ 
ing  points  in  the  wind  jet.  We  do  not  know  the  relevant 
special  functions  for  four  nearby  turning  points,  but 
their  computation  is  likely  to  be  too  time-consuming 
anyway  for  our  purposes. 

So  instead  we  model  four  turning  points  using  the 
interpolated  ray  solution  for  two  turning  points.  Three 
examples,  each  with  a  different  horizontal  wavelength 
A,  are  shown  in  Fig.  4.  The  left  panels  show  vv,  and  the 
right  panels  show  the  corresponding  nr  profile.  For 
comparison,  m2  is  also  plotted  without  wind  curvature 
terms.  With  wind  curvature,  the  mz  curve  has  two 
minima  and  more  closely  straddles  the  nr  =  0  axis  in  the 
region  of  the  turning  points. 

The  initial  condition  for  these  calculations  is  an 
upgoing  plane  wave  of  unit  amplitude  at  z  =  20  km.  In 
the  upper  two  rows  of  Fig.  4,  the  numerical  solution  and 
the  interpolated  ray  solution  agree  closely. 

The  agreement  is  not  always  so  good.  This  is  illus¬ 
trated  in  the  bottom  row  of  Fig.  4,  where  the  interpo¬ 
lated  ray  solution  below  the  wind  jet  has  an  amplitude 
that  is  about  half  that  of  the  numerical  solution,  al¬ 
though  the  phases  agree  well.  There  are  only  two 
turning  points  in  this  case,  with  nr  not  quite  crossing 
over  to  positive  values  near  z  =  10  km.  The  reason  that 
the  two-turning-point  theory  is  not  accurate  for  this 
two-turning-point  case  is  that  it  assumes  a  locally  qua¬ 
dratic  shape  for  m 2  near  the  turning  points  rather  than 
the  quartic  shape  realized  here.  On  the  other  hand,  we 
should  not  have  expected  the  two-turning-point  theory 
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Fig.  3.  Examples  of  the  interpolated  ray  solution  for  vv,  without  wind  curvature  effects.  Turning-point  heights  are 
indicated  by  horizontal  dotted  lines,  (top)  Same  cases  as  in  the  respective  panels  of  Fig.  2.  (bottom  left)  This 
illustrates  a  relatively  large  error  in  the  interpolated  ray  solution,  showing  the  imaginary  part  of  w  for  a  horizontal 
wavelength  is  14  km.  (bottom  right)  The  slowly  varying  parameter  e  of  (28)  for  the  solution  in  the  bottom  left  panel. 


to  be  accurate  for  four-turning-point  cases  either,  con¬ 
trary  to  what  we  find  in  the  upper  two  rows  of  Fig.  4. 
More  work  is  needed  to  determine  the  factors  that  af¬ 
fect  the  accuracy  of  two-turning-point  theory  in  these 
cases. 

In  the  rest  of  this  paper,  we  always  use  the  two- 
turning-point  interpolated  ray  solution,  despite  its  poor 
accuracy  in  some  cases. 


5.  Other  effects  needed  for  mountain-wave 
simulations 

a.  Wave  reflection  from  the  ground 

Mountain  waves  that  reflect  from  a  turning  point 
return  to  the  ground,  where  they  reflect  upward  again. 
We  have  not  yet  taken  into  account  the  contribution 
from  these  ground-reflected  waves.  A  way  to  do  so  is 
given  in  Broutman  et  al.  (2006)  and  is  summarized 
next. 


Each  reflection  from  the  ground  generates  a  new  in¬ 
cident  wave,  with  a  phase  shift  and  a  change  in  ampli¬ 
tude  relative  to  the  original  incident  wave.  The  phase 
shift  is  due  to  the  total  phase  change  in  propagating 
from  the  ground  to  the  lower  turning  point  and  then 
back  to  the  ground  again,  including  the  phase  shift  of 
—  tt/2  at  the  turning  point  and  the  phase  shift  of  ir  for 
reflection  at  the  ground.  The  total  phase  shift  is  then 

<t>  =  2  I  m  dz  +  tt/2  —  2cr,  (30) 

Jo 

where  a  is  defined  in  (11).  As  noted  below  (11),  u  is 
small  and  only  has  a  slight  effect  (at  most  a  few  percent) 
on  the  results  below. 

For  two  turning  points,  the  change  in  amplitude  is 
given  by 

R=\wr\!\wi\  =  eQ/  [e2Q  +  1)1/2.  (31) 
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Fig.  4.  Examples  of  the  interpolated  ray  solution  for  w,  with  wind  curvature  effects:  (left)  solutions  for  the  real  part 
of  tv;  (right)  profiles  of  m2  (normalized  by  k2)  corresponding  to  the  left  panels.  Horizontal  dotted  lines  indicate 
turning-point  heights.  Included  for  comparison  is  the  m2  profile  without  wind  curvature  (thin  line).  The  horizontal 
wavelength  A  in  each  calculation  is  indicated  in  the  m 2  plot. 


This  reflection  coefficient  is  the  ratio  of  the  ray  ampli¬ 
tudes  for  the  reflected  and  incident  wave  solutions  of 
(9)-(10).  As  the  distance  between  the  two  turning  points 
widens,  Q  increases  and  R  approaches  unity. 

To  incorporate  nr  reflections  from  the  ground,  we 
take  the  solutions  for  w  obtained  above  (ray,  uniform, 
numerical)  and  multiply  each  by 


nr 

S=J2Rnein‘s>.  (32) 

n  =  0 

The  number  of  ground  reflections  nr  depends  on  the 
time  at  which  we  evaluate  the  solution  and  can  be 
computed  from  the  group  velocity  for  each  ray.  The 
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WRF  Fourier  with  wind  curvature  Fourier  without  wind  curvature 
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Fig.  5.  The  solution  w(x,  z)  in  m  s_1  at  t  =  4  h  for  the  sech-squared  wind  profile:  (left)  WRF;  (middle)  interpolated  ray  with  wind 
curvature;  (right)  interpolated  ray  without  wind  curvature.  Florizontal  dotted  lines  indicate  the  turning-point  locations  for  the  resonant 
modes  given  in  (38)-(41).  The  mountain  profile  (36)  is  also  shown,  amplified  by  a  factor  of  4  for  visibility. 


time  that  it  takes  a  wave  group  to  propagate  from  the 
ground  to  the  lower  turning  point  and  back  to  the 
ground  is 

T  =  2£  ctd^  (33) 

where  cg  is  the  vertical  group  velocity.  At  time  f,  the 
number  of  ground  reflections  can  be  approximated  as 
tIT  rounded  to  the  nearest  integer.  This  approach  was 
tested  for  an  idealized  case  by  Broutman  et  al.  (2006) 
and  was  used  in  the  Jan  Mayen  simulation  of  EB06. 

b.  Waves  with  one  or  no  turning  points 

The  solution  for  a  Fourier  component  with  only  one 
turning  point  involves  the  Airy  function.  It  is  given  in 
the  appendix  of  EB06  and  will  not  be  repeated  here. 
The  solution  for  a  Fourier  component  without  turning 
points  is 

w  =  —  ia)h(mo/m)1/2eiq,  (34) 

where  q  =  fn  nidz  and  m0  is  m  at  z  =  0.  This  satisfies 
the  lower  boundary  condition  (3). 

c.  Anelastic  scaling 

So  far  we  have  made  the  Boussinesq  approximation. 
From  here  on,  we  include  anelastic  effects  with  the 
usual  density  scaling  of  the  previously  derived  expres¬ 
sions: 

v^(po/p)1/2h>,  (35) 

where  p  (z)  =  p0  exp(— z  /  ! hi)  is  the  mean  density  and  p0 
refers  to  its  value  at  z  =  0.  We  use  the  density  scale 
height  I  hi  =  7.5  km.  Anelastic  corrections  to  the  dis¬ 
persion  relation  (2),  involving  Hd ,  were  found  to  be 
unimportant  in  the  following  examples.  The  results  in 


this  paper  are  based  on  the  Boussinesq  dispersion  re¬ 
lation  (2). 

6.  Mountain  waves  in  a  sech-squared  wind 

We  first  consider  a  two-dimensional  example,  with 
wind  shear  given  by  the  profile  (26)  and  a  constant  N  = 
0.015  s-1.  The  topography  is 

h  =  h0/{l+x2/L2),  (36) 

with  ho  =  100  m  and  L  =  2.5  km. 

The  Fourier  solution  is  computed  from  (8)-(10)  and 
(29)  for  the  trapped  waves  and  from  (34)  for  the  prop¬ 
agating  waves.  The  grid  has  512  points  in  x,  with  a  grid 
spacing  of  1  km,  and  101  points  in  z  with  a  spacing  of  200 
m.  The  solution  is  calculated  at  f  =  4  h.  The  time  affects 
the  number  of  ground  reflections  for  the  trapped  waves; 
see  nr  in  (32). 

We  also  show  a  numerical  solution  using  the  Weather 
Research  and  Forecasting  (WRF)  model,  a  mesoscale 
model  described  in  Skamarock  et  al.  (2005).  The  grid  for 
this  computation  has  480  grid  points  in  x  with  a  grid 
spacing  of  250  m,  and  301  grid  points  in  z  with  a  spacing 
100  m.  A  sponge  layer  starts  at  a  height  of  20  km. 

Figure  5  shows  the  (left)  WRF  solution  and  (middle), 
(right)  Fourier  solution  with  and  without  wind  curva¬ 
ture  effects,  respectively.  Wind  curvature  is  included  or 
omitted  by  setting  a  =  1  or  a  =  0,  respectively,  in  (2). 
The  horizontal  lines  are  turning-point  heights  for  the 
resonant  modes  and  are  discussed  below. 

When  wind  curvature  is  included,  there  are,  out  of  a 
total  of  256  horizontal  wavenumbers,  40  for  which  the 
waves  are  vertically  propagating  without  turning  points. 
There  are  80  horizontal  wavenumbers  for  which  the 
waves  have  turning  points,  and  of  these,  5  have  four 
turning  points  and  75  have  two  turning  points.  (The 
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other  horizontal  wavenumbers  correspond  to  evanes¬ 
cent  waves  and  are  ignored.)  The  waves  with  four 
turning  points  have  horizontal  wavelengths  in  the  range 
of  about  11.4-12.5  km.  The  waves  with  two  turning 
points  have  shorter  horizontal  wavelengths.  When  wind 
curvature  is  omitted,  these  numbers  change  to  30 
propagating  and  90  trapped  horizontal  wavenumbers. 
All  trapped  waves  then  have  two  turning  points  and 
horizontal  wavelengths  less  than  16.5  km. 

Of  the  two  Fourier  solutions,  the  one  with  wind  cur¬ 
vature  (Fig.  5,  middle)  is  in  better  agreement  with  WRF. 
The  correlation  coefficient  for  the  WRF  and  Fourier 
solutions  increases  from  about  0.57  without  wind  cur¬ 
vature  to  about  0.86  with  wind  curvature. 

The  differences  between  the  Fourier  solutions  with 
and  without  wind  curvature  can  be  explained  by  con¬ 
sidering  the  resonant  modes.  These  are  defined  here  as 
those  waves  that  have  a  phase  change  <F  in  (30)  that  is  a 
multiple  of  2tt.  For  such  Fourier  components  there  is 
perfect  constructive  interference  between  the  associ¬ 
ated  rays  [i.e.  those  leaving  the  ground  for  the  first  time 
and  those  repeating  their  paths  (in  Fourier  space)  after 
each  ground  reflection]. 

For  the  Fourier  solution  with  wind  curvature,  there 
are  two  resonant  modes  resolved  by  the  present  grid. 
Their  horizontal  wavelengths  Ares  and  the  correspond¬ 
ing  (lower)  turning-point  heights  zres  were  computed 
numerically  to  be  (in  kilometers) 


~  5.1, 7.5, 

(37) 

~  3. 8, 5.8,  and 

(38) 

=  7,  7. 

(39) 

Here,  nr  refers  to  the  number  of  ground  reflections  at 
t  =  4  h,  as  used  in  (32). 

For  the  Fourier  solution  without  wind  curvature, 
there  are  three  resonant  modes: 


~  5.0,  7.1,  15.9, 

(40) 

~  3.9,  6.0,  9.2,  and 

(41) 

II 

o 

p\ 

u> 

(42) 

The  heights  zres  are  indicated  by  the  horizontal  dashed 
lines  in  Fig.  5.  For  the  WRF  solution,  the  heights  in  (38) 
are  plotted. 

The  presence  of  a  third  resonant  mode  with  Ares  ~  15.9 
km  in  the  case  without  wind  curvature  seems  to  explain 
the  differences  in  the  two  Fourier  solutions  of  Fig.  5. 
At  heights  below  6  km,  the  Fourier  solution  without 
wind  curvature  shows  longer  horizontal  scales  than  the 


Fourier  solution  with  wind  curvature,  consistent  with  its 
longer  resonant  mode  and  with  the  interference  of 
three,  rather  than  two,  resonant  modes.  At  heights 
above  6  km,  the  Fourier  solution  without  wind  curva¬ 
ture  has  larger  wave  amplitudes  because  the  third 
resonant  mode  has  a  turning-point  height  near  9-km 
altitude.  Because  this  turning  point  is  close  to  the  peak 
of  the  wind  jet,  there  is  significant  transmission  of  this 
resonant  mode  through  the  wind  jet. 

The  ratio  of  transmitted  to  incident  wave  amplitude  is 
T  =  [1  +  exp(2<2)]~1/z.  This  is  derived  from  the  ray 
solution  (8),  (9),  with  Q  given  by  (7).  Thus,  small  values 
of  Q  mean  significant  transmission.  The  values  of  Q 
corresponding  to  the  three  resonant  modes  in  (40)  are, 
respectively, 

Q=  f*2  \m\  dz  —  11,  5.2,  0.14.  (43) 

Jzi 

Of  these,  only  the  15. 9-km  mode  has  a  significant  value 
of  T  ~  0.65.  At  z  =  20  km  in  the  solution  without  wind 
curvature  (right  panel  of  Fig.  5)  the  first  four  downwind 
crests  are  separated  by  about  15  km,  close  to  the  reso¬ 
nant  value. 

For  the  solution  with  wind  curvature  (middle  panel 
of  Fig.  5),  the  downwind  crests  at  z  =  20  km  are  sep¬ 
arated  by  the  slightly  smaller  distance  of  about  12-14 
km.  This  may  be  related  to  the  fact  that  adding  wind 
curvature  shortens  the  maximum  horizontal  wave¬ 
length  for  the  trapped  solutions.  Without  wind  curva¬ 
ture,  this  maximum  wavelength  is  about  16.5  km.  With 
wind  curvature,  the  maximum  wavelength  falls  to  12.5 
km. 

Above  the  wind  jet,  and  downwind  from  the  moun¬ 
tain,  there  is  also  a  contribution  to  the  wavefield  from 
waves  that  just  miss  having  turning  points.  These  waves 
are  advected  downwind  as  they  refract  to  a  small 
wavenumber  near  the  tip  of  the  wind  jet.  They  have  less 
than  half  the  amplitude  of  the  contribution  from  waves 
that  tunnel  through  the  wind  jet.  We  can  see  this  by 
separately  plotting  the  wavefield  due  to  the  propagating 
waves  and  due  to  the  trapped  waves.  This  is  shown  in 
Fig.  6. 

7.  Mountain  waves  over  Jan  Mayen 

We  next  consider  a  three-dimensional  case  of  moun¬ 
tain  waves  generated  by  the  island  of  Jan  Mayen  in  the 
far  North  Atlantic.  The  problem  was  studied  by  EB06 
using  a  Fourier  solution  without  two-turning-point  ef¬ 
fects.  As  in  EB06,  we  use  radiosonde  winds  and  strati¬ 
fication  for  25  January  2000  and  realistic  topography. 
There  is  a  wind  jet  centered  at  about  10-km  altitude. 
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Fig.  6.  The  vertical  velocity  w(x,  z)  in  m  s-1  for  (left)  vertically  propagating  waves  and  (right)  vertically  trapped 
waves,  for  the  Fourier  solution  without  wind  curvature.  Horizontal  dashed  lines  indicate  turning-point  heights  for  the 
resonant  modes. 


The  zonal  wind  component  is  plotted  in  the  right  panel 
of  Fig.  1.  For  more  details,  see  EB06. 

We  computed  the  Fourier  solution  on  a  grid  with  1024 
points  in  x,  512  points  in  y,  and  a  grid  spacing  of  1  km  in 
each  direction.  For  this  grid  and  the  radiosonde  wind 
and  stratification  profiles,  there  are  about  37  000  Fou¬ 
rier  components  that  correspond  to  stationary  mountain 
waves.  Of  these,  about  23  000  have  at  least  one  turning 
point,  and  of  these  about  21  000  have  two  or  more 
turning  points.  These  numbers  vary,  depending  on  how 
much  smoothing  is  applied  to  the  background  wind  and 
stratification  profiles  and  on  how  the  background  wind 
and  stratification  is  approximated  at  heights  below  the 
mountain  top. 

We  limit  the  number  of  two-turning-point  solutions 
we  need  to  calculate  by  considering  the  integral  Q  in  (7). 
As  evident  from  the  discussion  below  (43),  large  values 
of  Q  correspond  to  waves  that  are  mostly  reflected  from 
the  wind  jet,  with  little  transmission  through  the  wind 
jet.  Thus  we  compute  two-turning-point  solutions  only 
for  waves  with  Q  <  3.  Horizontal  wavenumbers  with 
higher  values  of  Q  are  treated  here  as  one-turning-point 
problems,  using  an  Airy  function  solution  up  to  a  height 
midway  between  the  two  turning  points.  Above  that 
height,  these  wavenumbers  are  ignored  because  the 
exponential  decay  between  the  turning  points  is  large  by 
that  point.  This  is  the  procedure  used  in  EB06  for  all 
horizontal  wavenumbers  with  two  or  more  turning 
points. 

There  are  about  6000  horizontal  wavenumbers  with 
Q  <  3.  These  are  treated  with  the  interpolated  ray 
solution  (29).  The  three-dimensional  solution  is  calcu¬ 
lated  from  the  ground  to  a  height  of  20  km  at  1-km 
intervals.  Only  one  vertical  cross  section  is  shown  here. 


The  total  computation  time  on  a  2-GHz  processor  is 
about  3  min,  including  just  under  1  min  for  the  two- 
turning-point  calculation. 

The  Jan  Mayen  solutions  of  EB06  are  shown  in  the 
top  row  of  Fig.  7.  These  were  computed  by  (upper  left) 
the  mesoscale  model  WRF  and  (upper  right)  the  Fou¬ 
rier  method  without  two-turning-point  effects.  The  is¬ 
land  of  Jan  Mayen  is  centered  at  the  origin  of  the 
horizontal  axis. 

The  Fourier  solution  based  on  the  interpolated  ray 
solution  is  shown  in  the  lower  left  panel  of  Fig.  7.  Here 
we  see  that  there  is  much  better  agreement  with  the 
WRF  results  when  nearby  turning-point  effects  are  in¬ 
cluded  in  the  Fourier  method,  both  in  the  amplitude  of 
the  transmitted  waves  above  the  wind  jet  and  in  the 
associated  reduction  of  the  wave  amplitude  with  dis¬ 
tance  below  the  wind  jet.  Wind  curvature,  included  in 
the  calculation  for  the  lower  left  panel,  also  has  a  no¬ 
ticeable  effect.  The  Fourier  solution  without  wind  cur¬ 
vature  is  shown  in  the  lower  right  panel  of  Fig.  7.  Note 
that  the  dominant  horizontal  wavelength  of  the  trapped 
waves  below  about  10-km  height  is  somewhat  shorter 
than  in  the  WRF  result  or  in  the  Fourier  solution  with 
wind  curvature  (lower  left  panel). 

The  difference  between  the  calculations  with  and 
without  wind  curvature  in  Fig.  7  is  mainly  a  small  change 
in  the  dominant  wavelength  of  the  downwind  wavefield. 
Downwind  of  the  mountain  (at,  say,  z  =  5  km),  there  are 
about  four  and  one-half  wave  crests  in  the  calculation 
with  wind  curvature  (lower  left  panel)  and  five  crests  in 
the  calculation  without  wind  curvature  (lower  right 
panel).  In  three  dimensions,  there  is  a  continuum  of 
resonant  wavenumbers  in  k ,  l  space,  as  shown  in  Fig.  10 
of  EB06.  Thus,  small  changes  in  the  dispersion  relation 
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Fig.  7.  Vertical  cross  sections  of  the  vertical  velocity  w(x,  z)  in  m  s_1  at  t  =  4  h  for  mountain  waves  generated  by  the 
island  of  Jan  Mayen,  (top)  The  WRF  solution  and  Fourier  solution  of  EB06.  (bottom)  The  Fourier  solutions  that 
include  two  turning-point  effects,  both  (lower  left)  with  and  (lower  right)  without  wind  curvature.  The  horizontal 
distance  on  the  axis  refers  to  a  direction  that  is  36°  north  of  eastward  (see  Fig.  8  of  EB06). 


due  to  wind  curvature  lead,  in  general,  to  small  changes 
in  the  resonant  horizontal  wavenumber. 

In  two  dimensions,  as  in  our  earlier  example  of  Fig. 
5,  the  resonant  wavenumbers  are  discrete.  If  one 
of  those  resonant  wavenumbers  has  a  turning  point 
near  the  crest  of  the  wind  jet,  a  small  change  in  the 
dispersion  relation  can  eliminate  that  resonance.  The 
corresponding  rays  then  change  from  vertically  trap¬ 
ped  to  vertically  propagating  and  no  longer  contribute 
as  strongly  to  the  downwind  wavefield  below  the  wind 
jet. 

Thus,  in  our  examples,  we  find  a  difference  between 
the  effects  of  wind  curvature  in  two  and  three  dimen¬ 
sions.  In  two  dimensions,  where  the  resonant  wave- 
numbers  are  discrete,  the  effects  of  wind  curvature 
can  change  the  number  of  resonant  horizontal  wave- 
numbers.  In  three  dimensions,  where  the  resonant 
wavenumbers  are  continuous,  the  inclusion  of  wind 


curvature  slightly  changes  the  wavenumber  at  reso¬ 
nance. 

8.  Discussion 

We  dealt  with  the  breakdown  of  ray  theory  near 
closely  spaced  turning  points  in  a  wind  jet  by  linearly 
interpolating  the  ray  solution  across  the  turning-point 
region.  This  involved  ray  solutions  from  the  two-turn- 
ing-point  theory  of  Froman  and  Froman  (1970),  which 
were  used  here  even  when  there  were  four  turning 
points  in  the  wind  jet  region. 

For  cases  with  two  turning  points,  this  procedure 
worked  well,  as  illustrated  in  Fig.  3  for  a  sech-squared 
wind  profile.  The  success  of  the  interpolation  is  due  in 
part  to  having  a  reliable  criterion,  (27)  or  (28),  to 
identify  where  the  ray  solution  breaks  down  and  should 
be  replaced  by  the  interpolated  values. 
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For  cases  with  four  turning  points,  which  arise  when 
wind  curvature  terms  are  included  in  the  dispersion 
relation  (2),  this  procedure  produced  accurate  results 
for  some  Fourier  components.  This  is  illustrated  in  the 
top  two  rows  of  Fig.  4.  But  for  other  Fourier  compo¬ 
nents,  the  error  was  significant,  as  illustrated  in  the 
bottom  row  of  Fig.  4.  More  study  is  needed  to  under¬ 
stand  which  cases  work,  and  why  they  should  work  at 
all.  The  ray  solution  for  four  turning  points  derived  by 
Froman  and  Dammert  (1970)  may  be  helpful  in  this 
regard,  but  it  was  not  implemented  here. 

Another  area  for  further  study  is  waves  that  just  miss 
having  turning  points  in  the  wind  jet.  We  approximated 
these  Fourier  components  as  purely  upward-propagat¬ 
ing  waves  [see  (34)]  when  actually  they  would  generate 
significant  partial  reflection  from  the  wind  jet.  Thus,  the 
wave  amplitudes  for  these  Fourier  components  were 
probably  overestimated  above  the  wind  jet  and  under¬ 
estimated  below  the  wind  jet.  Partial  reflection  occurs  to 
some  extent  for  all  wave  propagation  in  a  variable  me¬ 
dium.  For  waves  that  just  miss  turning  points,  the 
turning  points  are  positioned  near  the  real  axis  in  the 
complex  z  plane,  and  the  strength  of  the  partial  reflec¬ 
tion  is  a  function  of  the  distance  from  these  complex 
turning  points  to  the  real  axis.  The  theory  is  given  in 
Froman  and  Froman  (1970). 

Despite  the  shortcomings  in  our  approach,  it  did  lead 
to  an  improved  comparison  with  the  results  of  a  meso- 
scale  model.  The  effects  of  nearby  turning  points  were 
important  in  both  examples,  as  were  the  effects  of  wind 
curvature  in  the  wind  jet. 

Other  methods  for  treating  nearby  turning  points 
may  turn  out  to  be  more  efficient.  Our  aim  here  was  to 
show  what  is  required  to  implement  ray  solutions  for 
nearby  turning  points.  We  have  not  experimented  at  all 
with  solutions  to  (1)  based  on  a  layered  approximation 
for  the  background  (e.g..  Smith  et  al.  2002;  Sutherland 
and  Yewchuk  2004;  Brown  and  Sutherland  2007),  and 
we  have  not  experimented  much  with  the  numerical 
integration  of  (1).  In  Nault  and  Sutherland  (2007),  a 
numerical  integration  of  an  equation  similar  to  (1) 
took  about  1  s  of  computation  time  per  Fourier  com¬ 
ponent  on  a  standard  processor,  and  accordingly  the 
integration  of  a  300  X  300  grid  of  Fourier  components 
took  about  1  day  to  complete.  However,  if  a  portion  of 
the  Fourier  spectrum  is  selected  such  that  the  same 
integration  step  size  is  appropriate  for  all  of  the  se¬ 
lected  Fourier  components,  the  integration  could  be 
vectorized  with  a  great  reduction  in  overall  computing 
time.  NaSu  have  begun  to  experiment  with  a  vector- 
izable  Heun’s  method  for  a  problem  with  nearby 
points.  Our  limited  tests  with  this  method  also  suggest 
that  it  might  be  effective.  More  testing  is  required,  and 


the  present  solutions  should  be  useful  for  comparison 
and  interpretation. 
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